#! /usr/bin/python

#
# Count numebr of reads mapping to regions hit by blast lincRNA vs pig genome
# BED file produced by 20110617_lincrna.sql
#

import sys
import os
import subprocess

os.chdir('/exports/work/vet_roslin_nextgen/dario/miscellanea/20110620_lincrna_pig')

bams=['/exports/work/vet_roslin_nextgen/dario/tophat/output/20110420_rnaseq_am_bmdm/am/ctrl/accepted_hits.bam',
      '/exports/work/vet_roslin_nextgen/dario/tophat/output/20110420_rnaseq_am_bmdm/am/lps/accepted_hits.bam',
      '/exports/work/vet_roslin_nextgen/dario/tophat/output/20110420_rnaseq_am_bmdm/bmdm/ctrl/accepted_hits.bam',
      '/exports/work/vet_roslin_nextgen/dario/tophat/output/20110420_rnaseq_am_bmdm/bmdm/lps/accepted_hits.bam']
ref_file= 'lincrna.blastn.bed'
outbeds= ['lincrna_am_ctrl.bed', 'lincrna_am_lps.bed', 'lincrna_bmdm_ctrl.bed', 'lincrna_bmdm_lps.bed']

for i in range(0, len(bams)):
    cmd= "coverageBed -abam %s -b %s > %s" %(bams[i], ref_file, outbeds[i])
    print(cmd)
    p= subprocess.Popen(cmd, shell=True, stdout= subprocess.PIPE, stderr= subprocess.PIPE)
    outerr= p.stderr.read()
    stdout= p.stdout.readlines()
    print(stdout)
    print(outerr)
